Improved Method for Detecting Local Discontinuities in CMB data by Finite Differencing 
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An unexpected distribution of temperatures in the CMB could be a sign of new physics. In partic- 
ular, the existence of cosmic defects could be indicated by temperature discontinuities via the Kaiser- 
Stebbins effect. In this paper, we show how performing finite differences on a CMB map, with the 
noise regularized in harmonic space, may expose such discontinuities, and we report the results of 
this process on the 7-year Wilkinson Microwave Anisotropy Probe data. 
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The Kaiser-Stebbins effect[l], a cosmological mani- 
festation of a discontinuous temperature gradient, is a 
well-known result of the presence of defects in the CMB 
such as cosmic strings[2]. The magnitude of the gradient 
relates to the string tension G\i by 
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where a s contains the dynamical and observational in- 
formation of the string such as the Lorentz factor and 
the orientation of the string segment; we will assume 
a s = 1. Current limits are of order G/i < 10 _6 [3]. 
One method of detecting such a discontinuity is by 
looking for characteristic signals in the derivatives of 
the CMB temperature field. Real-space derivative op- 
erators produce significant small-scale noise, limiting 
their utility for cosmology. This is particularly relevant 
when searching for highly localized sources such as lo- 
cal defects. The noise properties of harmonic and other 
transform-based methods, such as wavelets, tend to be 
well behaved up to the Nyquist limit but these trans- 
forms act to smooth out discontinuities. However, it is 
the different error characteristics of these two types of 
approaches (with respect to each other) in the presence 
of discontinuous signals that suggest one might calcu- 
late the difference of the derivative maps created by each 
method in order to isolate discontinuous signals. 
When calculating derivatives over a pixelated sampling 
grid, a finite difference scheme [4] is both popular and 
effective. It can be shown that a given finite difference 
scheme is related to an underlying interpolating poly- 
nomial [5]; this being the case, the scheme will perform 
poorly when in the presence of functions which are not 
well modeled by such a polynomial, due to Runges phe- 
nomenon [6]. This is true for discontinuous signals, with 
the error in the resultant derivative scaling with both the 
order of the differencing scheme and the severity (the 
step size) of the discontinuity. These errors can then be 
used as a crude detector of discontinuities. 
In the following text, a method to detect discontinu- 
ities and the results of an application to the Wilkinson 
Microwave Anisotropy Probe (WMAP) 7-year internal 
linear combination (ILC) CMB temperature anisotropy 
map are reported. The ILC map is known to suffer from 



systematic noise at multipoles I > 100, which we fur- 
ther confirm with the detection of a relic of the galaxy 
boundary. This is performed using the MASQU soft- 
ware package [5], which was initially designed for clean 
masked CMB polarization mode separation [7] on the 
popular HEALPIX [8] spherical grid, utilizing real-space 
differential operators which act on the basis-dependent 
orthogonal Stokes linear polarization parameters Q and 
U to produce the scalar and pseudoscalar e and b fields, 
real-space analogues of the E and B modes. In the past, 
the signal-to-noise ratio in the Stokes parameters has 
been too low for such a method; projects such as the 
Planck Surveyor [9] will go some way toward rectifying 
that. In the present case, derivative operations will be 
performed on the temperature anisotropy map instead. 
The derivatives of a field with respect to coordinate ba- 
sis n at each point on a pixelated grid can be computed 
by the finite difference approximation [4] 
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where the are numerical weights attached to a 
sample of surrounding pixels (the pixel stencil). The 
weights are calculated by inverting a matrix V at each 
pixel, whose elements depend on the positions of the 
stencil pixels. In the following, nth-order calculations 
O n refer to calculations using a square (n + l) 2 pixel 
stencil. Figure 1 shows how the results vary with both 
step size and stencil size, for the E/B mode-separation. 
Note that the error increases both in magnitude and in 
range for larger stencils in the presence of a disconti- 
nuity, whereas it would decrease for a function well- 
modeled by a polynomial. 

The ILC map[10] is an N S id e = 512 resolution tem- 
perature aniostropy map, created by the WMAP team 
by taking weighted combinations of the WMAP band- 
limited data (bands K,Ka,Q,V and W), with special at- 
tention paid to the masked region corresponding to the 
galactic plane. For our analysis of the ILC map, it is use- 
ful to calculate the Laplacian of the field, which corre- 
sponds to a simple power multiplier 
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FIG. 1: The effect of step size and stencil size on our deriva- 
tives. We calculate the pseudoscalar V 4 b on a function 
Q{6, 4>) = sin 9 cos 0, U = 0, with a cut-off to Q{9, 0) = at 
the equator (top-left diagram), where 6 is related to the po- 
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larization B-mode by b(h) 
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(n). Since 



there is no U signal, the b signal should be zero everywhere. 
Clockwise from top-right diagram: b fields for O2 , O4 and O g 
calculations. The map resolution here is N ai d e = 32, while the 
scale is logarithmic. 



when applied to the harmonic space reconstruction of a 
spherical function 
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calculated explicitly as 
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where the Yj TO are the ordinary spherical harmonics. 
From here on in, we refer to the method of taking a 
HEALPIX map, generating the a; TO s via the iterative 
HEALPIX map2alm technique[ll], multiplying these 
coefficients by the factor in Eq. (3), and creating a 
Laplacian map by the HEALPIX alm2map technique 
as merely the spectral method. The spectral method 
will also expose any conspicuous boundaries, and suf- 
fers from Gibbs' phenomenon[12,13], similar to that of 
Runge. Further, the iterative HEALPIX method used to 
recreate the ai m s is of reasonable but nevertheless lim- 
ited accuracy. These errors will also propagate through 
to the Laplacian map. Since the interpolating polyno- 
mial and the spherical harmonics are not identical, a dis- 
continuity should be enhanced by calculating the map 
of differences between the finite difference Laplacian 
and the spectral Laplacian. This is then particularly use- 
ful when the discontinuity is subtle. As a toy model for 
the calculation, we show the results for a CMB simula- 
tion injected with a diamond-shaped discontinuity (i.e. 
by linearly adding ST/T = 100/iK to the map at the dis- 
continuity region, see Fig. 2). 



FIG. 2: Top row, left-to-right: CMB simulation (T rms « 
90/iK)+diamond discontinuity, Laplacian map. The boundary 
of the discontinuity can be spotted by eye. Bottom row, left-to- 
right: difference maps for O2 and O4 MASQU calculations; the 
discontinuity boundary is clearly accentuated in these images. 
The by-eye detection limits for this model in the ILC map are 
ST/T ~ 40^K (0 2 ) and ~ 15^K (O4), which for strings would 
correspond to tensions of G[i ~ 1.5 x 10~ 6 and ~6x 1CT 7 



To estimate this calculation, it is noted that the differ- 
encing error for derivatives of a function f(0,(f>) at pixel 
i is approximately 
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where (n,p) are the orders of the derivatives in coordi 
nate bases (6, </>), yielding for the full Laplacian 
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Meanwhile, the function reconstruction error from 
the Gibbs' phenomena receives contributions from two 
places; First, the Fourier sum in the limit of I — > 00 will 
reconstruct the signal everywhere except at the discon- 
tinuity. Secondly, truncation creates overshoot (ringing) 
in the vicinity of the discontinuity. These errors will be 
of the form 
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where g is a truncated approximation to the original 
function /, cut off at l„ 
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FIG. 3: Thick line: power spectrum for the O2 absolute differ- 
ence map calculated over a map synthesized from a Gaussian 
random sample of et; m s provided by the (a = lfiK) HEALPIX 
random number generator facility. The dashed and dotted 
lines are the spectra for the spectrally-calculated and O2 finite- 
differenced Laplacian maps respectively. The discrepancy be- 
tween the dashed and dotted lines decreases with an increased 
number of sampling points in the finite-differencing calcula- 
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tions. The power spectrum C\ Bpoc F D of the map dif- 
ferences has a white spectrum on small scales due to the strong 
correlation of noise between the maps and a boost on large 
scales due to the HEALPIX geometry in the polar cap (see [5]). 

An analysis of this type on the sphere has been carried 
out for longitudinal ringing via the Gegenbauer poly- 
nomials in [14]. Finally, there is an error contribution 
<5/hpix which accounts for any small remaining errors 
accumulated from the HEALPIX routines. In calculat- 
ing the Laplacian of the field, the full harmonic function 
errors 5f spcc will have additional scaling as ~ I 2 . Our 
method thus calculates 

<5(V 2 /)., « 5(V 2 f) diS ,i - V 2 (<5/ spcCil ). (11) 

at each pixel position i. The noise characteristics (Fig. 3) 
of the absolute difference map accounts for why a dis- 
continuity can be seen by eye in the toy model absolute 
difference map; the ~ I 2 scaling in either Laplacian map 
is approximately canceled, leaving white noise on small 
scales. 

The results for the software performance on the ILC 
data are shown in Fig. 4. Taking the absolute difference 
between the finite-differenced ILC map and spectrally 
calculated ILC map reveals the galaxy mask boundary 
very clearly, regardless of the smoothing used in the ILC 
map construction[10]. 

The ILC map is not directly useful for cosmological 
parameter estimation; instead, it is used as an approx- 
imately noise-free CMB map in order to isolate galac- 
tic foregrounds and generate a zeroth-order noise tem- 
plate [15] which is then supplemented by a maximal 
entropy method. In this sense the impact of removing 



FIG. 4: Top row: ILC log Laplacian maps from the (left dia- 
gram) map2alm method and (right diagram) finite differenc- 
ing at O2 ■ Bottom row: (left diagram) the difference map be- 
tween the above. The galaxy mask is quite apparent, as are a 
couple of signals which are likely galactic foreground sources. 
The triangular parts in the polar cap are relics of the calcula- 
tion technique and HEALPIX geometry. The last plot is the 
mask map, with ILC technique cuttings shown. 

the boundary discontinuity region for analyses is un- 
likely to seriously effect the cosmological implications 
of the WMAP results. Instead of implementing the com- 
plicated WMAP-type analysis to reconstruct the angu- 
lar power spectrum, the boundary is isolated in a sep- 
arate mask. One can then estimate the contribution of 
the masked region to the ILC power spectrum by using 
the relation between the pseudo-C; spectrum C\ and the 
theoretical power spectrum given below: 

Ci = ^2k w C v (12) 

w 

via a mode-coupling matrix Kw . More typically, this 
technique is used in MASTER-type calculations[16] (a 
Monte Carlo-type power spectrum estimator, such as 
SPICE[17]) since the equation above properly recovers 
only the Monte Carlo average of the spectral signal. This 
is performed for a range of boundary thicknesses (2-,4- 
and 6-pixel mask thicknesses), since the smoothing ker- 
nel corresponds to 1.5dcg[18] which is not rigorously fol- 
lowed in the code used to create the boundary mask. 

It can be seen in Fig. 6 that the impact of the bound- 
ary region is small; however, as cosmological experi- 
ments probe the microwave sky with higher precision 
and theory discrimination becomes more sensitive to 
subtle differences in predicted power spectra, it may be 
necessary to ensure that excess power from the mask 
boundary does not bias results. Secondly, given recent 
claims of detections of non-Gaussianity[19,20] it would 
certainly be necessary to account for the boundary re- 
gion in any calculation of the nonlinear term /iVL[21]. 
The effect of the boundary region on a number of the 
so-called anomalies [22] has also been calculated, with 
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FIG. 5: Top row: Boundary masks used, of varying thickness. 
Bottom row: alteration in the quadrupole (left diagram) and 
octupole (right diagram) power from the smaller boundary 
mask. The alignment of the quadrupole and octupole is es- 
sentially unchanged. 
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FIG. 6: Differences between the boundary masked and un- 
masked ILC power spectra — the lines correspond to the thick- 
est mask (thick line), thinnest mask (dotted line) and interme- 
diate mask (dashed line). The dip in power at I ~ 180 cor- 
responds roughly to the Gaussian smoothing scale across the 
galaxy boundary in the construction of the ILC map, so the 
difference between ILC and mode-uncoupled masked spectra 
at the mask scale is small. 



no significant divergence from their manifestation when 
the boundary region is included (Fig. 5). 

One can continue the search for anomalies by exclud- 
ing the galaxy region precalculation and generating a 
new Laplacian-difference map. For firm statistical anal- 
ysis (eschewing any by-eye method), the needlet ap- 
proach^] is particularly suitable to decomposing a dis- 
continuity into power due to its localization properties, 
while the Canny algorithm[24] is suited to edge detec- 
tion. Further, it should be noted that there may be cor- 
relations at the discontinuity boundary with the Stokes' 
Q and U parameters, since cosmic strings have a vector 
S-mode signal[25] in real space analogous to the Kaiser- 



Stebbins effect in the temperature case. These issues will 
be discussed in a future paper. 

In conclusion, the method proposed finds no immediate 
evidence for anomalies such as strings, with the caveat 
that only lower-order calculations have been performed, 
coupled with by-eye detection. This finding may change 
with higher-order calculations and a more rigorous ex- 
amination using the needlet algorithm. Also, since the 
signal in a pixel is an average over signals within the 
pixel area, it is possible that the results of the Planck Sur- 
veyor with an improved resolution of iV^rfe = 2048, or 
more local balloon surveys such as the E and B Exper- 
iment [26] and their successors, will test these results. 
The method did encounter a success: the detection of 
the smoothed galaxy boundary in the map. It was found 
that the effect of excluding the boundary from the map 
was negligible, but that its mere detection might war- 
rant added care in future surveys when performing map 
synthesis. 
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